/*******************************************************************************

	SSC0143 - Programacao Concorrente
	Prof. Julio Cesar Estrela
	
	T1 - Calculo de PI com 10000000 de casas decimais
	
	Grupo: 
	
	Daniel de Oliveira Pereira		6427171
	Guilherme Takasawa Yagui		6426698
	Vinicius Marques Stocco			6427233
	
	Março de 2013
	
********************************************************************************	
				Algoritmo de Borwein - Versao Concorrente	
*******************************************************************************/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
#include <pthread.h>
#include <semaphore.h>

// numero de casas decimais corretas de PI
#define PREC_DGT	10000000

// numero de bits minimo para armazenar PI com o numero de casas decimais corretas 
#define PREC_BIT 	(int) ceil((log(10) / log(2)) * PREC_DGT)

// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado 
#define ITR 		(int) ceil(log(PREC_DGT) / log(4))

mpf_t ak, yk, ak1, yk1, bk;				// calculo de pi
mpf_t tmp1, tmp2, tmp3, tmp4;			// auxiliares
	
sem_t s_main, s_Y1, s_Y2, s_A1, s_A2; 	//semaforos
	
pthread_t t_y1, t_y2, t_a1, t_a2;		//threads

// calcula o numerador de y(k+1)
// y1 = (1 - (1 - y(k)^4) ^ (1/4))
void *Y1(void *trash)
{	
	while (1)
	{
		// auto-bloqueio
		sem_wait(&s_Y1);
	
		mpf_pow_ui(tmp1, yk, 4);
		mpf_ui_sub(tmp1, 1, tmp1);
		mpf_sqrt(tmp1, tmp1);
		mpf_sqrt(tmp1, tmp1);
		mpf_ui_sub(tmp1, 1, tmp1);
	
		// libera a thread principal
		sem_post(&s_main);
	}
}

// calcula o denominador de y(k+1)
// y2 = (1 + (1 - y(k)^4) ^ (1/4))
void *Y2(void *trash)
{
	while (1)
	{
		// auto-bloqueio
		sem_wait(&s_Y2);
	
		mpf_pow_ui(tmp2, yk, 4);
		mpf_ui_sub(tmp2, 1, tmp2);
		mpf_sqrt(tmp2, tmp2);
		mpf_sqrt(tmp2, tmp2);
		mpf_add_ui(tmp2, tmp2, 1);
	
		// libera a thread principal
		sem_post(&s_main);
	}
}

// calcula o termo 1 de a(k+1)
// a1 = ak * (1 + y(k+1)) ^ 4     
void *A1(void *trash)
{
	while (1)
	{
		// auto-bloqueio
		sem_wait(&s_A1);
	
		mpf_add_ui(tmp3, yk1, 1);
		mpf_pow_ui(tmp3, tmp3, 4);
		mpf_mul(tmp3, tmp3, ak);
	
		// libera a thread principal
		sem_post(&s_main);
	}
}

// calcula o termo 2 de a(k+1)
// a2 = y(k+1) * (1 + y(k+1) + y(k+1) ^ 2)
	void *A2(void *trash)
{
	while (1)
	{
		// auto-bloqueio
		sem_wait(&s_A2);
	
		mpf_pow_ui(tmp4, yk1, 2);
		mpf_add(tmp4, tmp4, yk1);
		mpf_add_ui(tmp4, tmp4, 1);
		mpf_mul(tmp4, tmp4, yk1);
		
		// libera a thread principal
		sem_post(&s_main);
	}
}

int main(char argc, char *argv[])
{
	FILE *f;
	
	clock_t t;		// tempo
	
	// define a precisão desejada para cada variavel
	mpf_set_default_prec(PREC_BIT);
	
	// marca  o clock inicial da execucao do algoritmo
	t = clock();
	
	f = fopen("borwein_paralelo.txt", "w");
	
	// aloca memoria para as variaveis utilizadas para calcular PI
	mpf_init(ak);
	mpf_init(yk);
	mpf_init(ak1);
	mpf_init(yk1);
	mpf_init(bk);
	mpf_init(tmp1);
	mpf_init(tmp2);
	mpf_init(tmp3);
	mpf_init(tmp4);
	
	// inicializa a(0)
	mpf_sqrt_ui(tmp1, 2);
	mpf_mul_ui(tmp1, tmp1, 4);
	mpf_ui_sub(ak, 6, tmp1);
	
	// inicializa y(0)
	mpf_sqrt_ui(tmp2, 2);
	mpf_sub_ui(yk, tmp2, 1);
	
	// inicializa os semaforos
	// libera Y1, Y2 e bloqueia A1, A2
	sem_init(&s_main, 0, 0);
	sem_init(&s_Y1, 0, 1);
	sem_init(&s_Y2, 0, 1);
	sem_init(&s_A1, 0, 0);
	sem_init(&s_A2, 0, 0);
	
	// cria as threads utilizadas para calcular y(k+1) e a(k+1)
	pthread_create(&t_y1, NULL, &Y1, NULL);
	pthread_create(&t_y2, NULL, &Y2, NULL);
	pthread_create(&t_a1, NULL, &A1, NULL);
	pthread_create(&t_a2, NULL, &A2, NULL);
	
	unsigned int i;
	unsigned int aux;
	
	for (i = 0; i < ITR; i++)
	{
		//aux 2 ^ (2k+3)
		aux = (unsigned int) pow(2.0, (2.0 * i + 3.0));
		
		// bloqueia a thread principal
		sem_wait(&s_main);
		sem_wait(&s_main);
		
		// calcula y(k+1)
		// y(k+1) = y1 / y2
		mpf_div(yk1, tmp1, tmp2);
		
		// libera as threads A1 e A2
		sem_post(&s_A1);
		sem_post(&s_A2);
		
		// bloqueia a thread principal
		sem_wait(&s_main);
		sem_wait(&s_main);
		
		// calcula a(k+1)
		// a(k+1) = a1 - aux * a2
		mpf_mul_ui(tmp4, tmp4, aux);
		mpf_sub(ak1, tmp3, tmp4);
		
		// calcula PI = b(k) = 1/a(k+1)
		mpf_ui_div(bk, 1, ak1);
		
		//gmp_printf("\n %.*Ff", PREC_DGT, bk);
		gmp_fprintf(f, "%.*Ff \n", PREC_DGT, bk);
		
		mpf_set(yk, yk1);
		mpf_set(ak, ak1);
		
		// libera as threads A1 e A2
		sem_post(&s_Y1);
		sem_post(&s_Y2);	
		
	}
	printf("\n Numero de iteracoes : %d \n", ITR);
	
	// cancela a execucao das threads
	pthread_cancel(t_y1);
	pthread_cancel(t_y2);
	pthread_cancel(t_a1);
	pthread_cancel(t_a2);
	
	// libera recursos alocados pelas threads
	pthread_join(t_y1, NULL);
	pthread_join(t_y2, NULL);
	pthread_join(t_a1, NULL);
	pthread_join(t_a2, NULL);
	
	// destroi os semaforos
	sem_destroy(&s_main);
	sem_destroy(&s_Y1);
	sem_destroy(&s_Y2);
	sem_destroy(&s_A1);
	sem_destroy(&s_A2);
	
	// libera a memoria alocada para calcular PI
	mpf_clear(ak);
	mpf_clear(yk);
	mpf_clear(ak1);
	mpf_clear(yk1);
	mpf_clear(bk);
	mpf_clear(tmp1);
	mpf_clear(tmp2);
	mpf_clear(tmp3);
	mpf_clear(tmp4);
	
	fclose(f);
	
	// calcula o numero clocks necessario para executar o algoritmo
	t = clock() - t; 
	
	printf("\n Tempo de execucao : %.3f s \n\n", (double) t / CLOCKS_PER_SEC);

	return 0;
}
